function result = gp(X, t, x_star, sigma=0.01, theta=0.1, l=10)
% there is not covariance function argument because octave does not support
% function pointers. 

% work with the mean of the target input (t)
mean_t = mean(t);
t -= mean_t;

% calculating big K
K = zeros(size(X, 1));
for x = 1:size(X, 1)    
    for y = 1:size(X, 1)
        K(x, y) = k(X(x, :)', X(y, :)', theta, l); 
    end
end

% calculating small k
k_star = zeros(size(X, 1), 1);
for i = 1:length(k_star)
    k_star(i) = k(X(i, :)', x_star', theta, l);
end

L = chol(K + sigma * eye(size(K)), 'lower');
alpha = L'\(L\t);
f = k_star'*alpha;
v = L\k_star;
sigma_star = k(x_star', x_star', theta, l) - v' * v;

n = size(L, 1);
LL = sum(L.*eye(size(L)));
log_p = -0.5*t'*alpha - sum(log(LL')) - (n/2 * log(2 * pi));

% dont forget to add the mean of t, which we removed at the beginning
result.F = f + mean_t;
result.SIGMA_STAR = sigma_star;
result.LOG_P = log_p;
